home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / linfit.pro < prev    next >
Text File  |  1997-07-08  |  6KB  |  150 lines

  1. ;$Id: linfit.pro,v 1.6 1997/01/15 03:11:50 ali Exp $
  2. ;
  3. ; Copyright (c) 1994-1997, Research Systems, Inc.  All rights reserved.
  4. ;       Unauthorized reproduction prohibited.
  5. ;+
  6. ; NAME:
  7. ;       LINFIT
  8. ;
  9. ; PURPOSE:
  10. ;       This function fits the paired data {X(i), Y(i)} to the linear model,
  11. ;       y = A + Bx, by minimizing the chi-square error statistic. The result
  12. ;       is a two-element vector containing the model parameters,[A,B].
  13. ;
  14. ; CATEGORY:
  15. ;       Statistics.
  16. ;
  17. ; CALLING SEQUENCE:
  18. ;       Result = LINFIT(X, Y)
  19. ;
  20. ; INPUTS:
  21. ;       X:    An n-element vector of type integer, float or double.
  22. ;
  23. ;       Y:    An n-element vector of type integer, float or double.
  24. ;
  25. ; KEYWORD PARAMETERS:
  26. ;   CHISQ:    Use this keyword to specify a named variable which returns the
  27. ;             chi-square error statistic as the sum of squared errors between
  28. ;             Y(i) and A + BX(i). If individual standard deviations are 
  29. ;             supplied, then the chi-square error statistic is computed as
  30. ;             the sum of squared errors divided by the standard deviations.
  31. ;
  32. ;  DOUBLE:    If set to a non-zero value, computations are done in double 
  33. ;             precision arithmetic.
  34. ;
  35. ;    PROB:    Use this keyword to specify a named variable which returns the
  36. ;             probability that the computed fit would have a value of CHISQR 
  37. ;             or greater. If PROB is greater than 0.1, the model parameters 
  38. ;             are "believable". If PROB is less than 0.1, the accuracy of the
  39. ;             model parameters is questionable.
  40. ;
  41. ;    SDEV:    An n-element vector of type integer, float or double that 
  42. ;             specifies the individual standard deviations for {X(i), Y(i)}.
  43. ;
  44. ;   SIGMA:    Use this keyword to specify a named variable which returns a 
  45. ;             two-element vector of probable uncertainties for the model par-
  46. ;             ameters, [SIG_A,SIG_B].
  47. ;
  48. ; EXAMPLE:
  49. ;       Define two n-element vectors of paired data.
  50. ;         x = [-3.20, 4.49, -1.66, 0.64, -2.43, -0.89, -0.12, 1.41, $
  51. ;               2.95, 2.18,  3.72, 5.26]
  52. ;         y = [-7.14, -1.30, -4.26, -1.90, -6.19, -3.98, -2.87, -1.66, $
  53. ;              -0.78, -2.61,  0.31,  1.74]
  54. ;       Define a vector of standard deviations with a constant value of 0.85
  55. ;         sdev = replicate(0.85, n_elements(x))
  56. ;       Compute the model parameters, A and B.
  57. ;         result = linfit(x, y, chisq = chisq, prob = prob, sdev = sdev)
  58. ;       The result should be the two-element vector:
  59. ;         [-3.44596, 0.867329]
  60. ;       The keyword parameters should be returned as:
  61. ;         chisq = 11.4998, prob = 0.319925
  62. ;
  63. ; REFERENCE:
  64. ;       Numerical Recipes, The Art of Scientific Computing (Second Edition)
  65. ;       Cambridge University Press
  66. ;       ISBN 0-521-43108-5
  67. ;
  68. ; MODIFICATION HISTORY:
  69. ;       Written by:  GGS, RSI, September 1994
  70. ;                    LINFIT is based on the routines: fit.c, gammq.c, gser.c,
  71. ;                    and gcf.c described in section 15.2 of Numerical Recipes,
  72. ;                    The Art of Scientific Computing (Second Edition), and is
  73. ;                    used by permission.
  74. ;         Modified:  SVP, RSI, June 1996
  75. ;             Changed SIG_AB to SIGMA to be consistant with the other
  76. ;             fitting functions. Changed CHISQR to CHISQ in the docs
  77. ;                    for the same reason. Note that the chisqr and the SIG_AB
  78. ;             keywords are left for backwards compatibility.
  79. ;         Modified:  GGS, RSI, October 1996
  80. ;                    Modified keyword checking and use of double precision.
  81. ;                    Added DOUBLE keyword.
  82. ;-
  83.  
  84. FUNCTION LinFit, x, y, chisqr = chisqr, Double = Double, prob = prob, $
  85.                        sdev = sdev, sig_ab = sig_ab, sigma = sigma
  86.  
  87.   ON_ERROR, 2
  88.  
  89.   TypeX = SIZE(X)
  90.   TypeY = SIZE(Y)
  91.   nX = TypeX[TypeX[0]+2]
  92.   nY = TypeY[TypeY[0]+2]
  93.  
  94.   if nX ne nY then $
  95.     MESSAGE, "X and Y must be vectors of equal length."
  96.  
  97.   ;If the DOUBLE keyword is not set then the internal precision and
  98.   ;result are identical to the type of input.
  99.   if N_ELEMENTS(Double) eq 0 then $
  100.     Double = (TypeX[TypeX[0]+1] eq 5 or TypeY[TypeY[0]+1] eq 5)
  101.  
  102.   nsdev = n_elements(sdev)
  103.  
  104.   if nsdev eq nX then begin ;Standard deviations are supplied.
  105.     wt = 1.0 / sdev^2
  106.     ss = TOTAL(wt, Double = Double)
  107.     sx = TOTAL(wt * x, Double = Double)
  108.     sy = TOTAL(wt * y, Double = Double)
  109.     t =  (x - sx/ss) / sdev
  110.     st2 = TOTAL(t^2, Double = Double)
  111.     b = TOTAL(t * y / sdev, Double = Double)
  112.   endif else if nsdev eq 0 then begin
  113.     ss = nX + 0.0
  114.     sx = TOTAL(x, Double = Double)
  115.     sy = TOTAL(y, Double = Double)
  116.     t = x - sx/ss
  117.     st2 = TOTAL(t^2, Double = Double)
  118.     b = TOTAL(t * y, Double = Double)
  119.   endif else $
  120.     MESSAGE, "sdev and x must be vectors of equal length."
  121.  
  122.   if Double eq 0 then begin 
  123.     st2 = FLOAT(st2) & b = FLOAT(b)
  124.   endif
  125.  
  126.   b = b / st2
  127.   a = (sy - sx * b) / ss
  128.   sdeva = SQRT((1.0 + sx * sx / (ss * st2)) / ss)
  129.   sdevb = SQRT(1.0 / st2)
  130.  
  131.   if nsdev ne 0 then begin
  132.     chisqr = TOTAL( ((y - a - b * x) / sdev)^2, Double = Double )
  133.     if Double eq 0 then chisqr = FLOAT(chisqr)
  134.     prob = 1 - IGAMMA(0.5*(nX-2), 0.5*chisqr)
  135.   endif else begin
  136.     chisqr = TOTAL( (y - a - b * x)^2, Double = Double )
  137.     if Double eq 0 then chisqr = FLOAT(chisqr)
  138.     prob = chisqr * 0 + 1 ;Make prob same type as chisqr.
  139.     sdevdat = SQRT(chisqr / (nX-2))
  140.     sdeva = sdeva * sdevdat
  141.     sdevb = sdevb * sdevdat
  142.   endelse
  143.  
  144.   sig_ab = [sdeva, sdevb]
  145.   sigma = sig_ab
  146.  
  147.   if Double eq 0 then RETURN, FLOAT([a, b]) else RETURN, [a, b]
  148.  
  149. END
  150.